NMDS on understory vegetation in Coastal Barrens group (MA & LI)

Libraries:

# libraries -------------------
library(ggplot2)
library(vegan)
library(ggvegan)
library(tidyverse)

Set working directory

knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')

Functions:

# functions -------------------
source("Functions/veganCovEllipse.R")

Source, tranform, select, and relativize MA & LI veg data (hidden)

Species found in at least 5% (11 plots) of plots in Coastal Barrens are: VAPA, GABA, QUIL, GAPR, VAAN, PTAQ, QUCO, PIRI, KAAN, MELI, RUSP, CAPE, QUAL, SMGL, & QUPR

Source and transform envr & categorical data; merge for predictor df to compare nmds against (hidden)

Run NMDS

# 3 dimensions
mali_3d <- metaMDS(rel_mali5,
                   distance = 'bray',
                   k = 3,
                   trymax = 20,
                   autotransform = FALSE)

## Run 0 stress 0.1272962 
## Run 1 stress 0.127296 
## ... New best solution
## ... Procrustes: rmse 0.0008255444  max resid 0.002723831 
## ... Similar to previous best
## Run 2 stress 0.1272962 
## ... Procrustes: rmse 0.0007830489  max resid 0.002595403 
## ... Similar to previous best
## Run 3 stress 0.1272961 
## ... Procrustes: rmse 0.0007857632  max resid 0.002554142 
## ... Similar to previous best
## Run 4 stress 0.1272959 
## ... New best solution
## ... Procrustes: rmse 0.0006474293  max resid 0.00211959 
## ... Similar to previous best
## Run 5 stress 0.1275742 
## ... Procrustes: rmse 0.03456454  max resid 0.1266969 
## Run 6 stress 0.1272961 
## ... Procrustes: rmse 0.0001239676  max resid 0.0004235327 
## ... Similar to previous best
## Run 7 stress 0.1466649 
## Run 8 stress 0.1582653 
## Run 9 stress 0.1275739 
## ... Procrustes: rmse 0.03522942  max resid 0.128812 
## Run 10 stress 0.1319007 
## Run 11 stress 0.1275737 
## ... Procrustes: rmse 0.0347356  max resid 0.1272635 
## Run 12 stress 0.1275743 
## ... Procrustes: rmse 0.03530544  max resid 0.1290684 
## Run 13 stress 0.127296 
## ... Procrustes: rmse 0.0001789619  max resid 0.000570497 
## ... Similar to previous best
## Run 14 stress 0.127296 
## ... Procrustes: rmse 9.379602e-05  max resid 0.0003106422 
## ... Similar to previous best
## Run 15 stress 0.1275734 
## ... Procrustes: rmse 0.03498567  max resid 0.1280777 
## Run 16 stress 0.1272959 
## ... New best solution
## ... Procrustes: rmse 7.379182e-05  max resid 0.0002473518 
## ... Similar to previous best
## Run 17 stress 0.127296 
## ... Procrustes: rmse 0.0001408249  max resid 0.0004555808 
## ... Similar to previous best
## Run 18 stress 0.1272962 
## ... Procrustes: rmse 0.0002469237  max resid 0.0008214091 
## ... Similar to previous best
## Run 19 stress 0.1275735 
## ... Procrustes: rmse 0.03493132  max resid 0.1278801 
## Run 20 stress 0.1275737 
## ... Procrustes: rmse 0.03467607  max resid 0.1270179 
## *** Best solution repeated 3 times

# converges, stress of 0.127 - 0.128

set.seed(2)
# 2 dimensions
mali_2d <- metaMDS(rel_mali5, 
                   distance = 'bray', 
                   k = 2,
                   try = 40,
                   trymax = 40,
                   autotransform = FALSE)

## Run 0 stress 0.2013098 
## Run 1 stress 0.2421891 
## Run 2 stress 0.2164813 
## Run 3 stress 0.2248403 
## Run 4 stress 0.2013098 
## ... New best solution
## ... Procrustes: rmse 7.336217e-05  max resid 0.0002404534 
## ... Similar to previous best
## Run 5 stress 0.2166909 
## Run 6 stress 0.2013098 
## ... Procrustes: rmse 2.20954e-05  max resid 7.187511e-05 
## ... Similar to previous best
## Run 7 stress 0.2299679 
## Run 8 stress 0.2334747 
## Run 9 stress 0.2067192 
## Run 10 stress 0.2063475 
## Run 11 stress 0.2248403 
## Run 12 stress 0.2052525 
## Run 13 stress 0.2013098 
## ... Procrustes: rmse 8.623799e-05  max resid 0.0002840637 
## ... Similar to previous best
## Run 14 stress 0.2248403 
## Run 15 stress 0.2013098 
## ... Procrustes: rmse 7.687422e-05  max resid 0.0002349284 
## ... Similar to previous best
## Run 16 stress 0.2385195 
## Run 17 stress 0.2378598 
## Run 18 stress 0.2013098 
## ... Procrustes: rmse 7.160904e-05  max resid 0.0002354649 
## ... Similar to previous best
## Run 19 stress 0.240332 
## Run 20 stress 0.271914 
## Run 21 stress 0.2063472 
## Run 22 stress 0.2013098 
## ... Procrustes: rmse 3.826719e-05  max resid 9.76465e-05 
## ... Similar to previous best
## Run 23 stress 0.2337596 
## Run 24 stress 0.2518641 
## Run 25 stress 0.2401184 
## Run 26 stress 0.2013098 
## ... Procrustes: rmse 6.309312e-06  max resid 1.438008e-05 
## ... Similar to previous best
## Run 27 stress 0.2431912 
## Run 28 stress 0.2013098 
## ... Procrustes: rmse 2.79146e-06  max resid 8.73469e-06 
## ... Similar to previous best
## Run 29 stress 0.2076391 
## Run 30 stress 0.2013099 
## ... Procrustes: rmse 0.0001601721  max resid 0.0005288505 
## ... Similar to previous best
## Run 31 stress 0.2013098 
## ... Procrustes: rmse 2.72951e-05  max resid 6.009303e-05 
## ... Similar to previous best
## Run 32 stress 0.2185542 
## Run 33 stress 0.2248403 
## Run 34 stress 0.2248403 
## Run 35 stress 0.239287 
## Run 36 stress 0.2248403 
## Run 37 stress 0.2076443 
## Run 38 stress 0.2394146 
## Run 39 stress 0.2013098 
## ... Procrustes: rmse 1.065455e-05  max resid 3.451836e-05 
## ... Similar to previous best
## Run 40 stress 0.2013098 
## ... Procrustes: rmse 2.107096e-05  max resid 6.965666e-05 
## ... Similar to previous best
## *** Best solution repeated 12 times

# stress still too high to represent in 2 dimensions, lowest is ~ 0.2

rm(mali_2d)

Make correlation data for bi-plots

ma.veg_enr <- envfit(mali_3d, ma.meta.data_veg)
ma.veg_spc.envr <- envfit(mali_3d, rel_mali5)

# site data -------------------
ma.site.scrs <- as.data.frame(scores(mali_3d, display = "sites"))

ma.site.scrs <- cbind(ma.site.scrs, Treat_Type = ma.meta.data_veg$Treat_Type)
ma.site.scrs <- cbind(ma.site.scrs, Region = ma.meta.data_veg$Region)

# species data   -------------------
ma.spp.scrs <- as.data.frame(scores(ma.veg_spc.envr, display = "vectors"))
ma.spp.scrs <- cbind(ma.spp.scrs, Species = rownames(ma.spp.scrs))
ma.spp.scrs <- cbind(ma.spp.scrs, pval = ma.veg_spc.envr$vectors$pvals)

sig.ma_spp.scrs <- subset(ma.spp.scrs, pval<=0.05)

# envr data -------------------
ma.envr.scrs <- as.data.frame(scores(ma.veg_enr, display = "vectors"))
ma.envr.scrs <- cbind(ma.envr.scrs, env.variables = rownames(ma.envr.scrs))
ma.envr.scrs <- cbind(ma.envr.scrs, pval = ma.veg_enr$vectors$pvals)
ma.sig_envr.scrs <- subset(ma.envr.scrs, pval<=0.05)

Need to think about how to represent 3D

Really need to figure out how to represent in 2D: figure out the % variation represented by each axis and use the 2 that represent the most variation (can also do multiple 2D graphs)

library(plotly)

fort.mali3 <- fortify(mali_3d) %>% 
  filter(score == "sites")



plot_ly(x = fort.mali3$NMDS1, y = fort.mali3$NMDS2, z= fort.mali3$NMDS3, type = "scatter3d", mode = "markers", color = ma.meta.data_veg$Treat_Type, symbols = ma.meta.data_veg$Region)

Run PerMANOVA to look at the impacts of region and treatment on community composition

summary(as.factor(ma.meta.data_veg$Treat_Type))

##  Control   FallRx  Harvest    MowRx SpringRx 
##        5        3        6        3        6

ma.tt <- adonis2(rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)

print(ma.tt)

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)
##            Df SumOfSqs      R2      F Pr(>F)    
## Treat_Type  4   1.9569 0.31903 2.1082  0.001 ***
## Residual   18   4.1769 0.68097                  
## Total      22   6.1339 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise follow up

library(pairwiseAdonis)

pairwise.adonis2(rel_mali5 ~ Treat_Type, data = ma.meta.data_veg)

## $parent_call
## [1] "rel_mali5 ~ Treat_Type , strata = Null , permutations 999"
## 
## $FallRx_vs_SpringRx
##            Df SumOfSqs      R2      F Pr(>F)
## Treat_Type  1  0.21549 0.10956 0.8613  0.613
## Residual    7  1.75138 0.89044              
## Total       8  1.96687 1.00000              
## 
## $FallRx_vs_MowRx
##            Df SumOfSqs      R2      F Pr(>F)
## Treat_Type  1  0.48859 0.46035 3.4122    0.1
## Residual    4  0.57276 0.53965              
## Total       5  1.06134 1.00000              
## 
## $FallRx_vs_Harvest
##            Df SumOfSqs      R2      F Pr(>F)
## Treat_Type  1  0.40379 0.20243 1.7767  0.109
## Residual    7  1.59094 0.79757              
## Total       8  1.99473 1.00000              
## 
## $FallRx_vs_Control
##            Df SumOfSqs      R2      F Pr(>F)
## Treat_Type  1  0.39033 0.23712 1.8649  0.118
## Residual    6  1.25581 0.76288              
## Total       7  1.64614 1.00000              
## 
## $SpringRx_vs_MowRx
##            Df SumOfSqs      R2      F Pr(>F)  
## Treat_Type  1  0.65697 0.28336 2.7678  0.015 *
## Residual    7  1.66151 0.71664                
## Total       8  2.31848 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $SpringRx_vs_Harvest
##            Df SumOfSqs      R2      F Pr(>F)  
## Treat_Type  1   0.4952 0.15598 1.8481  0.061 .
## Residual   10   2.6797 0.84402                
## Total      11   3.1749 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $SpringRx_vs_Control
##            Df SumOfSqs     R2      F Pr(>F)
## Treat_Type  1  0.16729 0.0666 0.6422  0.778
## Residual    9  2.34457 0.9334              
## Total      10  2.51186 1.0000              
## 
## $MowRx_vs_Harvest
##            Df SumOfSqs      R2      F Pr(>F)  
## Treat_Type  1  0.72529 0.32578 3.3823  0.014 *
## Residual    7  1.50107 0.67422                
## Total       8  2.22636 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $MowRx_vs_Control
##            Df SumOfSqs      R2      F Pr(>F)  
## Treat_Type  1  0.84147 0.41918 4.3303  0.014 *
## Residual    6  1.16594 0.58082                
## Total       7  2.00741 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $Harvest_vs_Control
##            Df SumOfSqs      R2      F Pr(>F)  
## Treat_Type  1  0.58418 0.21102 2.4072  0.023 *
## Residual    9  2.18412 0.78898                
## Total      10  2.76830 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## attr(,"class")
## [1] "pwadstrata" "list"

#significant: SpringRx vs. MowRx; SpringRx vs. Harvest (borderline at 0.051); MowRx vs. Harvest; MowRx vs Control; Harvest vs Control

Indicator species analysis

library(labdsv)

ma.meta.data_veg$Treat_Group <- NA

ma.meta.data_veg$Treat_Group <- ifelse( ma.meta.data_veg$Treat_Type == "FallRx", 1, ma.meta.data_veg$Treat_Group)

ma.meta.data_veg$Treat_Group <- ifelse( ma.meta.data_veg$Treat_Type == "MowRx", 2, ma.meta.data_veg$Treat_Group)

ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "SpringRx", 3, ma.meta.data_veg$Treat_Group)

ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "Harvest", 4, ma.meta.data_veg$Treat_Group)

ma.meta.data_veg$Treat_Group<- ifelse( ma.meta.data_veg$Treat_Type == "Control", 5, ma.meta.data_veg$Treat_Group)

mali_isa <- indval(mali5_avg, ma.meta.data_veg$Treat_Group)

summary(mali_isa)

##      cluster indicator_value probability
## QUPR       2          0.8549       0.002
## KAAN       2          0.7534       0.008
## RUSP       2          0.7097       0.005
## SMGL       4          0.6896       0.017
## 
## Sum of probabilities                 =  4.713 
## 
## Sum of Indicator Values              =  6.68 
## 
## Sum of Significant Indicator Values  =  3.01 
## 
## Number of Significant Indicators     =  4 
## 
## Significant Indicator Distribution
## 
## 2 4 
## 3 1

gr <- mali_isa$maxcls[mali_isa$pval <= 0.05]
iv <- mali_isa$indcls[mali_isa$pval <= 0.05]
pv <- mali_isa$pval[mali_isa$pval <= 0.05]
fr <- apply(mali5_avg > 0, 2, sum)[mali_isa$pval <= 0.05]
fridg <- data.frame(group=gr, indval=iv, pvalue=pv, freq=fr)
fridg <- fridg[order(fridg$group, -fridg$indval),]

print(fridg)

##      group    indval pvalue freq
## QUPR     2 0.8549332  0.002    6
## KAAN     2 0.7534400  0.008   10
## RUSP     2 0.7096551  0.005    7
## SMGL     4 0.6896051  0.017    6

# 1 is FallRx
# 2 is MowRx
# 3 is SpringRx
# 4 is Harvest
# 5 is Control

NMDS on understory vegetation in Coastal Barrens group (MA & LI)

Kathleen Stutzman

2024-02-26